'''
Created on Feb 2, 2010

@author: oabalbin
'''

import os
import sys
import glob
import shutil
import copy
import numpy as np
from optparse import OptionParser
from datetime import datetime
from collections import deque



# My libraries
import signatures.parsers.read_microarray_data as pm
import signatures.parsers.read_gene_lists as gpm
import signatures.db.tables as dt
import signatures.db.query as dq
import signatures.common.classes as cm #record, arrayMatrix
import signatures.preprocess.tools as tl #filter_constant_genes, get_variable_genes, 
import signatures.seedgenes.seedlist as sd #get_seed_genes, write_varin_file
import signatures.bfrm.paramfile as bp
import signatures.bfrm.bfrm_run as br
import signatures.bfrm.output_parser as bo
import signatures.plotting.plots as myplt


class pipeline():
    
    def __init__(self, dbhost, dbuser, dbpasswd, database):
        
        self.mt = dt.tables(dbhost, dbuser, dbpasswd, database)
        self.mq = dq.query(dbhost, dbuser, dbpasswd, database)
        #table information
        self.tableDescription = 'probeID VARCHAR(40), geneSym VARCHAR(40), sample VARCHAR(40), log_fold FLOAT(20,14), INDEX probe (probeID), INDEX gene (geneSym), INDEX sample (sample)'
        # table with in house prostate cohort
        #self.tableName = 'ERG_Prostate'
        
        # table with golub prostate cohort
        self.tableName = 'golub_prostate_cohort2'
        
        #self.tableName = 'NCI60_2'
        #self.tableName = 'NCI60'
        
        
        
        ######## Parameters ############
        # Parameters for pre-processing data:
        # percentage of low variable gene to be filter out of the analysis
        self.lvar=5 #15
        # percentage of low median gene to be filter out of the analysis
        self.lmed=5 #15
        
        
        # Parameters to run bfrm
        self.NLatentFactors = 2
        self.EvolIncludeVariableThreshold=0.75
        self.EvolIncludeFactorThreshold=0.75
        self.EvolMaximumFactors=0
        self.EvolMaximumVariables=0
        
        # Parameters to parse brfm output. 
        # inter= intercept (true or false), covar= number of covariates, nresfac=number of response factors, nlatfac= number of latent factors, 
        # ppith =  threshold. It determines a threshold for the probability of gene i to be included in factor j.
        self.inter, self.covar, self.nresfac, self.nlatfac, self.ppith = 0, 0 ,0, 0, 0.99 
        
        
        # sample classes
        
        self.spclass0=['RWPE_3_RWPE_7','RWPE_7_CPP','CPP_RWPE19','RWPE_3_CPP','RWPE_CPP_dup','CPP_RWPE_7',
                       'GEaN23','GEaN25','GEaN27','GEaN29','GEaN31','GEaN32','GEaN33',
                       
                       'GEaT32','GEaT38','GEaT40','GEaT42','GEaT43','GEaT45','GEaT47','GEaT49','GEaT53',
                       'GEaT54','GEaT55','GEaT56','GEaT57','GEaT58','GEaT61','GEaT65','GEaT66',
                       'GEaM16','GEaM20','GEaM26','GEaM33','GEaM36','GEaM39',
                       'GEaM40','GEaM43','GEaM44']
        

        self.spclass1=['PREC_AD_ERG1_LACZ','PREC_AD_ERG1_LACZ_dup','RWPE_L_ERG1_GUS','RWPE_L_ERG1_GUS_dup',
                       'PREC_AD_ETV1_LACZ','PREC_AD_ETV1_LACZ_dup','RWPE_L_ETV1_GUS','RWPE_L_ETV1_GUS_dup',
                       'PREC_AD_ERG1_LACZ_dyeswap','PREC_AD_ERG1_LACZ_dyeswap_dup','PREC_AD_ETV1_LACZ_dyeswap',
                       'PREC_AD_ETV1_LACZ_dyeswap_dup','RWPE_L_ERG1_GUS_dyeswap','RWPE_L_ERG1_GUS_dyes_dup',
                       'RWPE_L_ETV1_GUS_dyeswap','RWPE_L_ETV1_GUS_dyeswap_dup',
                       'GEaT20','GEaT52','GEaT62','GEaT64','GEaT76','GEaT77','GEaM15','GEaM21','GEaM22','GEaM23',
                       'GEaM28','GEaM29','GEaM34','GEaM37','GEaM38','GEaM41']
        
        self.spcltest=[]
        '''

        self.spcltest=['PREC_AD_ETV1_LACZ','PREC_AD_ETV1_LACZ_dyeswap','RWPE_L_ETV1_GUS_dyeswap','RWPE_L_ETV1_GUS','PREC_AD_ETV1_LACZ_dup','RWPE_L_ETV1_GUS_dup',
                     'PREC_AD_ETV1_LACZ_dyeswap_dup','RWPE_L_ETV1_GUS_dyeswap_dup','HPAEpiC_RAM_ERG_dup','HMEpiC_RAM_ERG_dup']
        
        self.spclass1=['PREC_AD_ERG1_LACZ','RWPE_L_ERG1_GUS','PREC_AD_ERG1_LACZ_dup','PREC_AD_ERG1_LACZ_dyeswap','RWPE_L_ERG1_GUS_dyes_dup',
                      'RWPE_L_ERG1_GUS_dyeswap', 'RWPE_L_ERG1_GUS_dup','PREC_AD_ERG1_LACZ_dyeswap_dup','PC3_ERG_dETS_1','PC3_ERG_dETS_3','PREC_RAM_ERG_dup','PREC_RAM_ERG']

        self.spcltest=['PREC_AD_ETV1_LACZ','PREC_AD_ETV1_LACZ_dyeswap','RWPE_L_ETV1_GUS_dyeswap','RWPE_L_ETV1_GUS','PREC_AD_ETV1_LACZ_dup','RWPE_L_ETV1_GUS_dup',
                     'PREC_AD_ETV1_LACZ_dyeswap_dup','RWPE_L_ETV1_GUS_dyeswap_dup','HPAEpiC_RAM_ERG_dup','HMEpiC_RAM_ERG_dup','HPAEpiC_RAM_ERG','HMEpiC_RAM_ERG']

        '''

        
    
    
    def parse_input_annotfiles(self, annotInputFile):
    
        expr_file, seed_file, outdirs, db_file, samples, bfrm, param, typefile={},{},{},{},{},{},{},'nofile'
    
        for line in annotInputFile:
            fields = line.strip('\n')
                 
            if fields[0][0] == '' or fields[0][0] == '#':
                continue
            
            if fields[0][0] == '@':
                typefile = fields
                continue
            
            fields = fields.split('=')
            dict = {'@Array expression':expr_file, '@Seed genes':seed_file,'@Database file':db_file,'@Samples':samples,'@Output Directories':outdirs,\
                    '@BFRM':bfrm,'@Parameters':param}[typefile]
            if typefile=='@Parameters':
                dict[fields[0]] = fields[1].split(',')
            else:
                dict[fields[0]] = fields[1]
            
            
        return expr_file, seed_file, db_file, samples, bfrm, outdirs, param
    

    def read_files_folder(self,folderpath,ext):
        ''' '''
        # Read files in folder
        myfiles=[]
        for infile in glob.glob( os.path.join(folderpath, '*'+ext) ):
            myfiles.append(infile)
            #print "current file is: " + infile
        return myfiles
    
    
    def process_microarray_files(self,exp_files,tp,soft=False, matfile=''):
        dbf=[]
        for f in exp_files:
            f1 = f.split('/')[-1]
            of=f+'_dbfile'
            if f1 == 'nci60_RNA_Affy_HGU95_AE_RMA_sp.txt':
            #if f == '/home/oabalbin/projects/networks/exp_data_house/prostate_cell_lines_ERG/raw_data/nci60_RNA_Affy_HGU95_AE_RMA_sp.txt':
                print 'reading nci60 samples ' + f
                tp.microarray_nci60_file(open(f),open(of,'w'))
            elif soft:
                print 'reading soft samples ' + f
                sp=f+'_sampleslist'
                gene=f+'_genelist'
                spstatus = f+'_sampleslist_status'
                sp_tr=f+'_sampleslist_training'
                sp_tt=f+'_sampleslist_test'
                trstatus = f+'_sampleslist_training_status'
                ttstatus = f+'_sampleslist_test_status'
                
                tp.microarray_exp_softfile(open(f),open(of,'w'),open(gene,'w'), open(sp,'w'), open(spstatus,'w'), open(sp_tr,'w'), \
                                           open(sp_tt,'w'), open(trstatus,'w'), open(ttstatus,'w') )
            elif matfile:
                print 'reading a GEO matrix file ' + f
                tp.microarray_exp_geoMatrixfile(open(f), open(matfile), open(of,'w'))
            else:
                print 'reading house samples ' + f
                tp.microarray_exp_file(open(f),open(of,'w'))
            dbf.append(of)
        
        return dbf
    
    def create_database(self,dbfiles,doit=False):
        """
        It creates a databases if doit = True
        """
        if doit:
            print 'Making tables'            
            for dbf in dbfiles:     
                if not(self.mt.table_existence(self.tableName)):
                    self.mt.create_table(self.tableName,self.tableDescription)
                    
                self.mt.load_data_table(self.tableName,dbf)
                
    def initialize_database(self, exp_files,tp, doit=False, soft=False, matfile=''):
        """
        Initialize and create the tables in the database
        """
        dbfiles = self.process_microarray_files(exp_files,tp, soft, matfile)
        self.create_database(dbfiles,doit)
    
    def check_create_dir( self, full_path_name ):
        """
        creates an output folder. Name folder is the date it was created
        returns the full path to the folder and the folder name only: outfolder, t.strftime
        """
        t = datetime.now()
        outfolder = full_path_name + t.strftime("%Y_%m_%d_%H_%M")
        if not os.path.isdir(outfolder):
            os.mkdir(outfolder)
        
        return outfolder, t.strftime("%Y_%m_%d_%H_%M")
    
    def move_bfrm_files(self, folderpath, tstampedfolder):
        """
        moves .txt files created for bfrm to the output time stamped directory
        """
        
        outtxtfiles = self.read_files_folder(folderpath,'.txt')
        for f in outtxtfiles:
            print f
            shutil.move(f,tstampedfolder)
    
    def check_bfrm_output(self, bfrmfolder,newFFile=[]):
        """
        Check if files  mA.txt, mF.txt, mPostPib.txt, exist
        """
        if newFFile:
            print newFFile
            bfrmof = deque(['mA.txt', 'mPostPib.txt', 'mF.txt','mPsi.txt','mVariablesIn.txt',newFFile])
        else:
            bfrmof = deque(['mA.txt', 'mPostPib.txt', 'mF.txt','mPsi.txt','mVariablesIn.txt'])
        
        files=[]   
        
        for f in bfrmof:
            fname = bfrmfolder + f
            print fname
            if os.path.isfile(fname):
                files.append(fname)
   
        return files
     
    def write_expmat_sdv(self, expMat, genelistind, outfile1, outfile2, spclass0, spclass1=[], spcltest=[]):
        """
        writes a expression matrix. The matrix is ordered according to the sample clasess: spclass0, spclass1, spcltest
        expMat= an expression Array object (common classes), genelistind =  the index in expMat of the wanted 
        genes.  string outfile1= writting file for the gene list, outfile2=writting file for the expression matrix
        """
        #self.mq.create_connector()
        self.mq.write_variable_genes_sdv(expMat, genelistind, outfile1, outfile2, spclass0, spclass1, spcltest)
    
    def write_arraytofile(self,outfile, H):
        """
        It writes the and numpy array to a file. 
        """
        for i in range(H.shape[0]):
            el = map(str,list(H[i,:]))
            outfile.write(",".join(el).replace(',','\t')+'\n')
            
    def write_bfrm_Hfile(self, outfile, H):
        """
        It is used to write the H array for bfrm including the design variables. 
        As: 03-11-2010. No control variables are included for running BFRM yet.
        """
        self.write_arraytofile(outfile, H)
    
    def write_bfrm_Ffile(self, outfile1, outfile2, H,samplelist):
        """
        write the F file.
        """
        outfile2.write(",".join(samplelist).replace(',','\t')+'\n')
        self.write_arraytofile(outfile1, H)

    def calc_corr_compfactor(self,pathactmat,compoundsmat,cmpsamples,samplelist,plotsfolder, corr_param, outfile20):
        """
        pathactmat= path activity matrix file name, compoundsmat = compounds dump file.
        """
        # load path activity matrix and compound sensitivity matrix. skip first row of intercepts
        # mat1 is the drug matrix. mat2 is the activity matrix
        # Samples names in the compound data set
        mat1_hd = tl.get_header_mat(open(cmpsamples))
        # load the compound matrix object and the path activity
        mat1 = tl.pickleload_matrix_data(open(compoundsmat))
        mat2 = np.loadtxt(pathactmat,dtype=float, delimiter='\t') #tl.get_matrix_data(open(pathactmat))
        
        #Filter the pathway activity matrix to obtain only the data belonging to the NCI60 or other samples.
        #nci60sp = np.array(expMat.sampleInd.values())
        #mat2= mat2[:,nci60sp]
        
        # Rearrange the compound matrix to match the column order of the
        # activity matrix. This is required in order to do the correlations
        odmat1 = tl.order_arrayMat_bysamples(mat1, samplelist, mat1_hd)  
        '''
        print mat2.shape
        print odmat1.expVal.shape
        '''
        
        outa = tl.calc_correlation_mats(odmat1.expVal, mat2, odmat1.geneSym, corr_param[0], float(corr_param[1]), float(corr_param[2]), outfile20 )
        # pickle the correlation matrix
        print 'Dumping the correlation array'
        outa.dump(plotsfolder+'pathcorr_dump.txt')
        # Plot the correlation matrices
        myplt.write_drugpath_mat(odmat1.expVal,mat2, outa, outdirs['outdir'],odmat1.geneSym)
        myplt.plot_drugvsact(odmat1.expVal,mat2, outa, plotsfolder)




        

##########
if __name__ == '__main__':
    
    optionparser = OptionParser("usage: %prog [options] ")
    optionparser.add_option("-f", "--annotFile", dest="annotFile",
                            help="annotation file for all files to use")
    optionparser.add_option("-d", action="store_true", dest="dbtrue",
                            help="create a database. Bool")
    optionparser.add_option("-t", action="store_true", dest="softtrue",
                            help="use a microarray soft file. Bool")    
    optionparser.add_option("-p", "--platformFile", dest="platformFile",
                            help="use a microarray soft file. Bool")


    # Running parameters
    optionparser.add_option("-u", action="store_true", dest="allgenes",
                            help="use genes in all uscs. Default use genes in universe file. Default universe. Bool")
    optionparser.add_option("-n", action="store_true", dest="norm",
                            help="normalize log2 and scale expression vals. Default only log2.  Bool")
    optionparser.add_option("-m", action="store_true", dest="normals",
                            help="standarize by normal expression and scale expression vals. Default only log2.  Bool")
    optionparser.add_option("-s", action="store_true", dest="seed",
                            help="use all genes for bfrm. Default seed genes. Bool")
    optionparser.add_option("-c", action="store_true", dest="drugcorr",
                            help="performs drug correlation. Default do not do it. Bool")
    optionparser.add_option("-a", "--nfactors", dest="nfactors",type="int",default=20,
                            help="Set the max number of factors for BFRM.")
    optionparser.add_option("-v", "--nvariables", dest="nvariables",type="int", default=1000,
                            help="Set the max number of variables for BFRM.")



    
    (options, args) = optionparser.parse_args()
       
    # start the pipeline. 
    # read input parameters
    pipe = pipeline("localhost", "oabalbin", "oscar", "signatures")
    # Set values for BFRM
    pipe.EvolMaximumFactors = options.nfactors
    pipe.EvolMaximumVariables = options.nvariables
    # Get the option from the input_param file
    expr_folder, seed_file, db_file, samples, bfrm, outdirs, run_param = pipe.parse_input_annotfiles(open(options.annotFile))
    
    # Create experiment folder tree and log file. 
    odir = cm.directoryTree()
    tstampedfolder, tstamp = odir.create_tstampdir(outdirs['outdir'])
    #' /bfrm, /drugsen, /drugsen/plots,/sdv' 
    dirs = odir.make_directory_tree()
    logfile = cm.logFile(tstampedfolder+'/'+tstamp+'_readme.txt') 
    
    ####### output files
    
    outfile1 = dirs[0]+tstamp+'_dataset.mat.txt'  # Expression values
    outfile2 = dirs[0]+tstamp+'_dataset.gen.txt'  # geneNames
    outfile7 = dirs[0]+tstamp+'_dataset.sample.txt'  # sampleNames in bfrm
    outfile3 = dirs[0]+tstamp+'_varin.txt'        # seed genes
    outfile4 = dirs[0]+tstamp+'_parameters.txt' #'outfolder+'/'+tstamp+'_parameters.txt'   # bfrm parameters file
    outfile9 = dirs[0]+tstamp+'_H.txt'          # H file with the design variables
    outfile5 = dirs[0]+tstamp+'_geneInModules.txt'

    
    outfile6 = dirs[3]+'/'+tstamp+'_mat.sdv.test.txt'   # sdv exp file gene
    outfile8 = dirs[3]+'/'+tstamp+'_gen.sdv.test.txt'   # sdv exp file exp values
    outfile20 = dirs[1]+tstamp+'permutationsFile.txt'          # It tells were is the permutation file used for the drug permutation analysis
    
    
    ######### Get raw data  ############
    # Get the database files
    #read folder for expression profiles
    exp_files = pipe.read_files_folder(expr_folder['expdir'],'.txt')

    tp = pm.parser()
    pipe.mt.create_connector()
    gtp = gpm.geneparser()
    ######## Create database
    if options.dbtrue:
        pipe.initialize_database(exp_files,tp, options.dbtrue, options.softtrue, options.platformFile)
    #sys.exit(0)
        
    # list of samples and all potential genes to be included in the study based 
    # on a pre-compiled gene list. 
    #Note: if an unbiased initial analysis wanted, then you need to run using the ucsc
    # list of all ucsc potential genes to include in study
    #samplelist = tp.list_of_names_in_line(open(samples['samples']))
    samplelist, bfrm_Hfile = tp.list_of_names_group(open(samples['samples']))
    
        
    ######### Preprocess to the expMatrix regarding the samples to be included to run BFRM
    
    # Order the samples according to a particular classification: e.i, ETS status
    # This function is particular to ETS annotation. 
    #If additional annotation wants to be included, 
    # The function and the assignment should be modified.
    if 'spannot' in samples.keys():
        spannot = tp.list_of_samples_status(open(samples['spannot']))  
        pipe.spclass0=spannot['ETS+']
        pipe.spclass1=spannot['ETS-']
        pipe.spclass2=spannot['None']
        sampleclasses = [pipe.spclass0, pipe.spclass1, pipe.spclass2]
        samplelist = tl.get_reordered_samplelist(samplelist,sampleclasses)
    
    #print samplelist
    ##### Gene lists
    
    if options.allgenes:
        extgenlist = gtp.list_of_genes_kgXref(open(seed_file['allgene']))  
    else:
        extgenlist = tp.list_of_names(open(seed_file['univfile']))
    
    # Get seed genes
    seedgenes = tp.list_of_names(open(seed_file['seedfile']))
    # Re-order extgenlist putting the seed genes at the begining of it .
    # It is required for BFRM, because it uses the first genes as seed for factors.
    extgenlist = tl.get_reordered_geneList(extgenlist,seedgenes)
    
    
    ########## Make queries ############
    # Do query to obtain the expression matrix
    # list of external genes to be used as the "universe of genes"
    # e.i: this could be the ERG interactome
    pipe.mq.create_connector()
    
    # Find a common set of genes for all samples, including the other samples
    if 'othersamples' in samples.keys():
        other_samples = tp.list_of_names_in_line(open(samples['othersamples']))
        allsamplelists = copy.copy(samplelist)
        allsamplelists.extend(other_samples)
        extgenlist = pipe.mq.get_commongen2samples(allsamplelists, extgenlist, pipe.tableName)

    ### Create the expMat array
    expMat_bfrm = pipe.mq.get_exp_mat_2(samplelist, extgenlist, pipe.tableName) 
    # Create sample index
    expMat_bfrm.create_sample_index()
    # Normalize the expression values data. log2 and oncomine scaling, or only log2 transform.
    # or it normalize using quantile normalization
    
    if options.norm:
        #expMat_bfrm.normalization()
        expMat_bfrm.quantile_norm()
        #myplt.density_plot(dirs[3]+'/'+tstamp+'_qnorm_plot.png',expMat_bfrm.expVal.filled(fill_value=0))
    else:
        expMat_bfrm.log2_norm()
        #myplt.density_plot(dirs[3]+'/'+tstamp+'_density_plot.png',expMat_bfrm.expVal.filled(fill_value=0))
    '''
    if options.normals:
        normal_samples=tp.list_of_names_in_line(open(samples['normals']))
        expMat_bfrm.normalizeexp_bynormals(normal_samples)
    '''
    
    ######### Pre-process Raw Data ###############
    # NOTE: This parameter pipe.lvar=25 can be modified to control variability filter
    # Get the index of genes with var >=lvar, and med >=lmed     
    #if dofilter:
    # Find variable genes and modify he expMat.
    genelist = tl.filter_out_constant_genes(expMat_bfrm, pipe.lvar, pipe.lmed, True)       
    # Get the index and name of the variable genes for a particular experiment. 
    #It returns a list of tuples with index and gene and a gen dictionary
    
    mgenelist, dictmgl =  tl.get_variable_genes(expMat_bfrm,genelist)    #np.arange(len(expMat_bfrm.genInd))  genelist
    expMat_bfrm = tl.get_filtered_expMat(expMat_bfrm,mgenelist,len(samplelist))
    
    # Index of genes belonging to the "universe" which also have var .=lvar and med >=lmedS
    expMat_bfrm.create_gen_index()
    genInd = expMat_bfrm.genInd.values()
    
    ########## Get seed genes
    #seedgenes = tp.list_of_names(open(seed_file['seedfile']))
    seeding, seednotfound = sd.get_gen_index(expMat_bfrm,seedgenes)
    
    
    ######### If nici60 samples provided or other additional samples provided.
    # create a bfrm object for them using the nci60samples and the final genes in the
    # expMat_bfrm which is going to be used to run BFRM    
    
    
    if 'othersamples' in samples.keys():
        outfile10 = dirs[0]+tstamp+'_mF_othersp.txt'
        outfile11 = dirs[0]+tstamp+'_mF_othersp_samples.txt'
        # Get the NCI60 sample list
        #other_samples = tp.list_of_names_in_line(open(samples['othersamples']))

        expMat_othersp = pipe.mq.get_exp_mat_2(other_samples, expMat_bfrm.geneSym, pipe.tableName)
        # Create sample index
        expMat_othersp.create_sample_index()
        # Normalize the expression values data. log2 and oncomine scaling, or only log2 transform.
        if options.norm:
            # This only applies to the nci60 cell lines, because they already quantile normalized
            expMat_othersp.log2_norm()
            #expMat_bfrm.quantile_norm()
            #expMat_othersp.normalization()
        else:
            expMat_othersp.log2_norm()
    
    
    ######## Write expression matrix for variable genes
    ######## write in varin.txt file the index of seed of genes. To use in BRFM
    ######### BFRM
    if options.seed:
        # Choose as seed all genes in the universe
        sd.write_varin_file(expMat_bfrm.genInd.values(),open(outfile3,'w'))
        pipe.mq.write_variable_genes(expMat_bfrm, expMat_bfrm.genInd.values(), open(outfile1,'w'), open(outfile2,'w'))
        pipe.write_bfrm_Hfile(open(outfile9,'w'),bfrm_Hfile)
    else:        
        # Choose as seed the genes given by the seedfile
        sd.write_varin_file(seeding.values(),open(outfile3,'w'))
        pipe.mq.write_variable_genes_nfact(expMat_bfrm, expMat_bfrm.genInd.values(), open(outfile1,'w'), open(outfile2,'w'), open(outfile7,'w'), False)
        pipe.write_bfrm_Hfile(open(outfile9,'w'),bfrm_Hfile)
        #pipe.mq.write_variable_genes_nfact_nrespon(expMat_bfrm, expMat_bfrm.genInd.values(), open(outfile1,'w'), open(outfile2,'w'), open(outfile7,'w'),
        #                                            pipe.spclass0, pipe.spclass1, pipe.spcltest)
    
    
    #########  BFRM #######
    ######### Create the BFRM parameters file
    # Call to parameters function
    '''
    #parmeter_val(self, NLatentFactors, NObservations, NVariables, DataFile, EvolVarIn, EvolVarInFile,
    #              EvolIncludeVariableThreshold=0.75, EvolIncludeFactorThreshold=0.75,
    #              EvolMaximumFactors=10, EvolMaximumVariables=150):
    '''
    NLatentFactors = pipe.NLatentFactors
    NObservations = len(samplelist)
    NVariables = len(expMat_bfrm.genInd) # keeping the analysis at the gene level. If the analysis made at probe level  this should change
    NDesignVariables = bfrm_Hfile.shape[1]
    DataFile = outfile1
    HFile = outfile9
    EvolVarIn = len(seeding)
    EvolVarInFile = outfile3
    EvolIncludeVariableThreshold= pipe.EvolIncludeVariableThreshold
    EvolIncludeFactorThreshold= pipe.EvolIncludeFactorThreshold
    EvolMaximumFactors= pipe.EvolMaximumFactors
    EvolMaximumVariables= pipe.EvolMaximumVariables
    
    
    ####### Write to logFile
    logfile.this.write('Input File '+ str(options.annotFile)+'\n')
    logfile.this.write('Database creation '+ str(options.dbtrue)+'\n')
    logfile.this.write('allgenes '+ str(options.allgenes)+'\n')
    logfile.this.write('norm '+ str(options.norm)+'\n')
    logfile.this.write('seed '+ str(options.seed)+'\n')
    logfile.this.write('The number of genes to study ' + str(len(extgenlist))+'\n')
    logfile.this.write('The Number of variable genes is '+str(len(genelist))+'\n')
    logfile.this.write('The number of seed genes is ' + str(len(seedgenes))+'\n')
    logfile.this.write('The number of variable seed genes is ' + str(len(seeding))+'\n')
    logfile.this.write('BFRM EvolIncludeVariableThreshold ='+ str(EvolIncludeVariableThreshold) +'\n')
    logfile.this.write('BFRM EvolIncludeFactorThreshold ='+ str(EvolIncludeFactorThreshold) +'\n')
    logfile.this.write('BFRM EvolMaximumFactors ='+ str(EvolMaximumFactors)+'\n')
    logfile.this.write('BFRM EvolMaximumVariables ='+ str(EvolMaximumVariables)+'\n')
    logfile.this.write('Samples in experiment'+'\n')
    logfile.this.write('Samples File '+ samples['samples']+'\n')
    logfile.write_list_names(samplelist)
    logfile.this.write('Samples Classes'+'\n')
    if 'spannot' in samples.keys():        
        logfile.write_list_names(spannot.keys())
        [logfile.write_list_names(spclass) for cl, spclass in spannot.iteritems()]  
        logfile.this.write('Genes in experiment'+'\n')
        logfile.write_list_names(seeding.keys())

    
    ######### RUN BFRM #######
    bfbp = bp.paramfile()
    
    bfbp.parmeter_val(NLatentFactors, NObservations, NVariables, NDesignVariables, DataFile, HFile, EvolVarIn, EvolVarInFile, 
                      EvolIncludeVariableThreshold, EvolIncludeFactorThreshold, EvolMaximumFactors, EvolMaximumVariables)
    
    bfbp.write_bfrmparam_file(open(outfile4,'w'))
    
    #sys.exit(0)
    ########## Create an bfrm_run object 
    bfrm = br.bfrm_run(outfile4,bfrm['bfrm'])
    bfrm.execute_bfrm()
    # need to change to the bfrm folder
    pipe.move_bfrm_files(os.getcwd(), dirs[0])
    
    ######### Parse the bfrm output to find the more significant genes in each module
    # Check if the files mPostPibFile, mAfile and mF.txt are in directory. Returns bfrmof = ['mA.txt', 'mPostPib.txt', 'mF.txt']
    
    bfof = pipe.check_bfrm_output(dirs[0])
    
    if len(bfof) < 4:
        print 'Not all BFRM output files were found'
        sys.exit(0)
    else:
        mA = np.loadtxt(bfof[0],dtype=float, delimiter='\t')
        mPostPib = np.loadtxt(bfof[1],dtype=float, delimiter='\t')
        mPsi = np.loadtxt(bfof[3],dtype=float, delimiter='\t')
        mVarIn = np.loadtxt(bfof[4],dtype=int, delimiter='\t')
        #print mVarIn
        #hvec = np.ones(mPsi.shape[0])
        nvar, nfactors = mA.shape[0], mA.shape[1]
    
    bfrmparser = bo.bfrm_parser(pipe.inter, pipe.covar, pipe.nresfac, pipe.nlatfac, pipe.ppith, nvar, nfactors)
    # it returns a dictionary of the final genes, and thir ind, present in all factors
    # Usage: self, outfile, mPostPibFile, mAfile, genelist, geneseed
    #geneinfactor = bfrmparser.get_gene_in_module(open(outfile5,'w'), bfof[1], bfof[0], expMat_bfrm.geneSym, seeding.keys())
    print expMat_bfrm.geneSym
    print mVarIn-1
    geneinfactor = bfrmparser.get_gene_in_module(open(outfile5,'w'), mPostPib, mA, expMat_bfrm.geneSym, seeding.keys(), mVarIn-1)
        
    logfile.this.write('########### BFRM Output '+ '\n')
    logfile.this.write('BFRM probability of inclusion (ppith) ='+ str(pipe.ppith)+'\n')
    logfile.this.write('BFRM number of factors (nfactors) ='+ str(nfactors)+'\n')
    logfile.this.write('BFRM number of variables (nvar) ='+ str(nvar)+'\n')
    

    ########## Pathway activity calculations. 
    ### Project the tumor factor into the other dataset:
    if 'othersamples' in samples.keys():
        # project the bfrm_ factor into the new data set.
        # factscore = bf.calc_factorscore_bfrm( mA[:,1:], mPostPib[:,1:],expVal,mPsi,hvec, B )
        
        # In the case in which there is not regression factors B is mean expression of each gene accros the samples.
        # 03-14-2010. It should be the mean expression of the genes in the data set used in BFRM.
        # However I going to use the mean expression in the nci60 data set instead. 
        Bmat = np.ma.mean(expMat_othersp.expVal.filled(fill_value=0),axis=1)
        othersamples_Fmatrix = bfrmparser.calc_factorscore_bfrm(mA[:,1:], mPostPib[:,1:],expMat_othersp.expVal[mVarIn-1,:].filled(fill_value=0),
                                                                 mPsi, np.ones(expMat_othersp.expVal.shape[1]), Bmat[mVarIn-1], 
                                                                 True)
        #othersamples_Fmatrix = bfrmparser.calc_factorscore_bfrm(mA[:,1:], mPostPib[:,1:],expMat_bfrm.expVal[mVarIn-1,:].filled(fill_value=0),
        #                                                         mPsi, np.ones(expMat_bfrm.expVal.shape[1]), mA[:,0], 
        #                                                         True)
        pipe.write_bfrm_Ffile(open(outfile10,'w'),open(outfile11,'w'),othersamples_Fmatrix,other_samples)
    
    
    '''
    # Use to evaluate correlation between the factor scores from BFRM and calculated from the posteriors using the calc_factorscore function
    mF = np.loadtxt(bfof[2],dtype=float, delimiter='\t')
    corr1,pval = tl.calc_corelation_arr(othersamples_Fmatrix[0,:], mF[1,:], False)
    myplt.plot_scatter2(othersamples_Fmatrix[0,:], mF[1,:], [corr1,pval], dirs[2]+'/')
    corr1,pval = tl.calc_corelation_arr(othersamples_Fmatrix[1,:], mF[2,:], False)
    myplt.plot_scatter2(othersamples_Fmatrix[1,:], mF[2,:], [corr1,pval], dirs[2]+'/')
    '''

    
    #### 04-06-2010 Note update with the new function for heatmaps.....
    
    ### Heat maps
    # vg_expMat = expMat.expVal[vargenes,:]
    # Get the expression matrix for the pathways activity heat maps
    pipe.write_expmat_sdv(expMat_bfrm, geneinfactor.values(), open(outfile6,'w'), open(outfile8,'w'), samplelist, spclass1=[], spcltest=[])
    # dirs[3]= folder directory for plot,expMat = exp matrix,row_names gene=names, col_names=sample names
    # The expMat.exVal is a masked array, convert to normal array with masked values = 0
    myplt.plot_heatmap(dirs[2]+"/heatmap_from_python.png", expMat_bfrm.expVal.filled(fill_value=0),expMat_bfrm.geneSym, samplelist)
    
    
    ########### CORRELATION between pathway activity and drug sensitivity
    # Check if the compound matrix is in a particular folder and run the analysis using the file.
    # if the file is not there run the file.
    # Check for the command line that tells to run the drug correlation
    
    # corr_param a list of 3 values that define the type of correlation performed
    corr_param=run_param['corr_param']
    
    logfile.this.write('########### Path Activity and GI50 correlation '+ '\n')
    logfile.this.write('The correlation parameters were =  ' + str(corr_param[0]) +' '+str(corr_param[1])+' '+str(corr_param[2]))
    logfile.this.close()
    
        
    if options.drugcorr:
        if ( os.path.isfile(db_file['dbdrug']) and os.path.isfile(db_file['cmpsamp']) and os.path.isfile(outfile10) ):
            # pathactmat,compoundsmat,cmpsamples,samplelist
            pipe.calc_corr_compfactor(outfile10,db_file['dbdrug'],db_file['cmpsamp'],other_samples, dirs[2]+'/',corr_param, outfile20)
    else:
        sys.exit(0)
        
        
    #logfile.this.close()
    sys.exit(0)
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    
    